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LOW-ORDER CLASSICAL RUNGE-KUTTA FORMULAS 
WITH STEPS IZE CONTROL 

AND THEIR APPLICATION TO SOME HEAT TRANSFER PROBLEMS 

INTRODUCTION 


1. In an earlier report [1], the author derived fifth- to eighth-order RUNGE- 
KUTTA formulas with stepsize control. In this paper, similar first- to 
fourth-order formulas are developed. 

2. Such low-order RUNGE-KUTTA formulas are of interest in some heat 
transfer problems. It is well-known that the parabolic partial differential 
equations of such problems can be reduced to ordinary differential equations. 
For instance, by a discretization of the space variable(s) of the problem, 
we obtain a system of ordinary differential equations with the time as the 
independent variable. Such a system can then be integrated by RUNGE- 
KUTTA methods. 

3. However, it is also well-known that the application of RUNGE-KUTTA 
methods to such problems is often very time-consuming. Higher-order 
RUNGE-KUTTA formulas do not offer advantages in this respect, since 
stability considerations, resulting from the exponential character of the 
solution, exclude an increase of the integration stepsize that would make 
such high-order formulas meaningful. Therefore, low-order RUNGE- 
KUTTA formulas (second- or third-order) can be expected to solve such 
problems more efficiently than any high-order formula. On the other hand, 
they are potentially more efficient than the standard difference formulas 
obtained by discretization of the space variable (s) as well as the time 
variable. 

4. For the efficiency of RUNGE-KUTTA formulas, it is essential that their 
truncation errors be as small as possible, since the permissable integra- 
tion stepsize is strongly dependent upon the magnitude of these errors. 
Therefore, we have tried to establish RUNGE-KUTTA formulas with small 
truncation errors. 

SECTION I. FOURTH-ORDER FORMULAS 

5. We consider the (vector) differential equation 

y 1 - f(x, y) , ^ 

and write for our RUNGE-KUTTA formula , 



f o = 


(K = i, 2, 3, 4,5) 


and 

4 

y = y + h Y. cf + 0(h 5 ) 

O K K 

K= 0 

5 

y = y +h £ cf + 0(h e ) 

0 K =0 K K 

with h as integration steps ize and (x , y ) as initial values. Equations 
(3) imply that we try to determine the coefficients , /3 , , c^ in 

such a way that the first formula (3) represents a fourth-order, and the 
second formula (3) , a fifth-order RUNGE-KUTTA formula. The difference 
y - y then represents an approximation for the leading (fifth- order) trun- 
cation error term of our fourth-order RUNGE-KUTTA formula and can be 
used easily for establishing a reliable stepsize control procedure for this 
formula. 

6 The coefficients a , 3 , c , and c have to satisfy certain equations of 

k kX k k 

condition that can be obtained by TAYLOR expansions. These equations of 
condition are well-known in the literature, see for example J.C. BUTCHER'S 
paper ( [2] , Table f) , or a paper by this author ( [ 1] , Table I) . For the 
convenience of the reader, we list these equations here for a fifth-order 
formula like the second equation (3) . 


f = 




Introducing the abbreviations: 


+ *K2 a * + * 


A „ p 

OL — P 

KK~i K ~ i K X 


(4) 


(k = 2, 3, 4, 5; X = 1,2,3) 

Table I contains the 17 equations of condition for our fifth-order formula. 
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TABLE I. EQUATIONS OF CONDITION FOR FIFTH-ORDER FORMULA 
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TABLE I. (Concluded) 


( V, 3) 
(V,2) 
(V,$) 
( V, 6) 

A 

(V,7) 
( V, 8) 
( V, 9) 




40 


- 0 


i 

6 


Z 

k=2 


A -o 

C P 
/C k3 


1 

120 


= 0 



1. 

2 


i s 

,y-0 


A 

c a 

K K 



30 


= 0 



A t^2 
C P 

K k! 


40 


0 



A 

C 

K 





J_ y 5 A 4 1 

24 ^ C /< a /< 120 

K=t 


0 


The Roman numerals in front of the equations in Table I indicate the order 
of the terms in the TAYLOR expansion. 

A similar table holds for a fourth-order formula such as the first formula 
( 3) . We obtain this table from Table I by omitting the fifth-order equations 
(V, 1) through (V, 9) and replacing, in the remaining equations, c bye 

K K 

and the upper limit 5 of the k-sums by 4. These remaining eight equations 
might be denoted by (I, 1) through (IV, 4) . 
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All eight equations of this new table for a fourth -order formula and the 17 
equations of Table I for a fifth -order formula have to be satisfied simultan- 
eously. 

7. For the following we assume 

cj = 0 , c* = 0 , a 4 = 1 ( 5 ) 

We further assume that a 2 , a 3 , a 5 are different from one another and from 
0 and 1. 

Equations (II, 1) , (HI, 2) , (IV, 4) then yield 


1 ~ 1 

12 o >2 (a? 3 - ot 2 ) ( i - oi 2 ) 


_i 2 QL \ ~ 1 

12 a 3 (a 2 " a 3) ( 1 " «3) 


1 6ao<^3 ~ 4( o?2 + Qf 3 ) + 3 

12 ( 1 - a 2 ) ( 1 ” cx 3 ) 


and equations (II, 1) , (III, 2) , (IV, 4) , (V, 9) 




> (6) 


/ 


A 

c 2 


1 1 Pomaces ~ 5(o?3+as) + 3 

60 a 2 (a 3 ~ a 2 ) ( i “ a 2 ) (“5 ~ ^ 2 ) 




a _ 1 _ 10 ceaQ's - 5(0*2 + Of 5 ) + 3 

° 3 “ 60 a 3 (a 2 - ce 3 ) (1 - a 3 ) (a 5 - a 3 ) 


A 

C4 = 


1 300*20*30*5 - 20 (a 2 a 3 + q;?q*r + Q^s) ± 15(^2 + o *3 + ^5) 

60 (a* 2 ~ 1) (a 3 - 1) (a* 5 “ 1) 


- 12 


> (?) 


A _ _1_ 1 0q ! 2 Q* 3 ~ 5(Q?2 + Qgs) + 3 

° 5 60 « 5 (q !2 “ a 5 ) ( 0*3 - o' 5 ) ( 1 - oi 3 ) 


) 
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8. Furthermore, we make the following assumptions that greatly reduce the 
number of equations of condition 


(Al) 


P 


21 


P31 - 


1 2 
2 «2 

1 2 
77^3 


P 41 - o 


P 51 = 77 


2 


(A2) 


P 22 “ 


32 


P 42 = 


52 


1 3 

3« 2 


3« 3 


i_ 

3 

i 3 

3 


(B) C 2 /3 2 i + C 3/ 3 31 + c 4/^41 - 0 

(B) C 2 P 2 I + c 3^ 31 + c 4/^41 + c 5^51 = 0 


(6) 


C2« 2 /3 2 i + c 3 cr 3 /3 3i + C4^4j + C5C2 5 /3 51 - 0 
From (Al) and (A2) the following identities result: 


A 


A 


(111,1) 

= (111,2) ; 

(IV, 2) 

= (IV, 4) ; 

(IV, 3) 

s 3 (IV, 4) 

1 

A 

A 

A 

A 

A 

= 3 ( IV , 4) ; 

} (8) 

(m, i) 

= (III, 2) ; 

(IV, 2) 

= (IV, 4) ; 

(iv, 3) 

A 

A 

A 

A 

, A 

A 

I 

( V, 6) 

= 4(V, 9) ; 

(V, 7) 

s 3(V,9) ; 

(V, 8) = 

6 ( V, 9) 


• using 

A 

(B) and (B) 

we find 

the following identities : 


(IV, 1) 

= (IV, 2) 






A 

A 

A 

A 

A 

A 

(9) 

(IV, 1) 

= (IV, 2) ; 

(V, 2) 

= ( V, 4) ; 

( V, 3) = 

3 ( V , 4) 



A 


and finally by also taking into account assumption ( C) 


( V, 5) = (V,6) 


( 10 ) 


AAA 

Therefore, equations (1,1) ; (1,1) ; (V, 1) ; (V, 4) are the only remaining 
equations of condition. 
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The first two equations determine the coefficients c and c that otherwise 

o o 

do not enter our equations v: condition. 

A A 

The remaining equations (V, 1) and (V, 4) have to be solved together 
with our assumptions (Al) , (A 2 ) , (B) , (B) , and (6) . 

9. From the first equations (Al) and (A 2 ) , the following relations are obtained: 



2 

= g«2 

(11) 

^21 

3 

= 1“ 2 

(12) 


The second equations (Al) and (A 2) yield 



The third equations (Al) and (A 2) , together with (B) s determine the 
coefficients j3 4l9 (S i2 $ i3 43 . We find the following expressions for these 



Elimination of c 5 from (B) and (C) leads to the following relation: 


c 2 (q! 5 - Q!2)/ 5 21 + c 3 («5 ■ 3) /3 31 + c 4 ( - l)/3 4i - 0 (15) 

A A A 

Introducing the above computed values for c 2 , c 3 , c 4 , jS 21 , /3 31 , and /3 41 into 
(15) leads to a relation between a 2 , <23, and a 5 „ 
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M q ?5 = N 


with 


(16) 


M - ( 6 a 2 Q !3 “ 4a 2 - 4ce 3 + 3) (30a 2 <^3 “ 30a 2 a 3 “ 10a 3 + 5 a 2 a 3 ) 

+ ( 6 ce 2 ce 3 - 6 o ' 2 «3 + 2a 3 -a 2 )(30a 2 ol \ ~ 20a 2 ce 3 - 20a 3 + 15a 3 ) 


?(i7) 


N = ( 6a 2 ce 3 “ 4a 2 - 4a 3 + 3) ( 16a 2 a | ~ 1 5a 2 cu 3 - 6a | + 3a 2^3) 

9 9 9 

+ ( 6a 2 a 3 “ 6a 2 O' 3 + 2a 3 — a 2 ) ( 20a 2 a 3 — 1 5 a 2 a 3 — 1 5 a 3 + 1 2a 3)/ 
It is easily verified that 


M = 0 (18) 

for any value of a 2 and a 3 . 

Because of ( 16) , only such values as a 2 and a 3 are possible that also lead 
to 


N = 0 


(19) 


Equation (19) represents a restrictive relation between a 2 and a 3 . This 
relation can be reduced to 


— . n ^ 

2 5a - 4a 2 + 1 


( 20 ) 


10. We use relation (20) to eliminate a 3 from the expressions for our RUNGE- 
KUTTA coefficients as obtained in Nos. 7, 8 , and 9. The elimination results 
in 


«i = 


a 3 = 


c 2 = 


2 

1 . “2 

2 5a 2 — 4a 2 + 1 

2 

1_ ' 1 5a 2 ~ 5a 2 + 1 

6 a 2 ( 1 - a 2 ) 10a| ~ Sa 2 + 1 


I (21) 


/ 
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2 1 ( 5 a\ - 4 a 2 + 1) 3 


c 3 = 

— • 

3 

C2 2 ( 5g! 2 " 2) 10a 2 

C 4 = 

i 

6 

10a 2 - 12a 2 + 3 
(1 - a 2 ) ( 5 «2 “ 1 2 ) 

/3 21 - 

3 

4 

a 2 

fi 31 = 

3 

15q; 2 - 12a 2 + 2 

16 

_ct 2 ’ (5 q! 2 - 4a 2 + i)° 

032 = 

- 

1 10(22 - 8(2 2 + 1 
8 a<l ( 5 a 2 ~ 4a 2 + 1) ' 

041 = 

3 

4 

1 

10ce 2 - 12 a 2 + 3 




( 21 ) 

( continued) 


1 1 - q!2 6 Oof ^ ~~ 78& 2 + 3 la 2 j 4 

^42 = “ 2 * (10a ^ 2 - 12q ! 2 + 3) ( 10(22 - 8 C 2 2 + 1) 

(1 - go) ( 5 a 2 - 2 ) (2g? ? . - 1) (5o| - 4 « 2 + l ) 2 
^43 2 q'2(10q'2 _ - 12^2 + 3) (10q!2 8 a 2 + 1) 


11. The weight factors £ 2 through c 5 of the fifth-order formula can now be 
expressed by a 2 and <a 5 : 


A 

C 2 - 


60 


i0ar(5al 5q/ 9+ D ■ ( 30a; jj - 29a 2 + 6) 
o;'{( 10a 2 - 8<a 2 + 1)(1 " <a 2 ) («5 “ «2) 


\ 


A 

^3 


A 

C4 


4 ( 5 a| - 4 q! 9 + l ) 4 [5Q!fi( 2 a 2 - 1) - ( 5c2 2 3)] 

15 * «|( 10^1 - 8 a 2 + i)(5<^ - 2 ) ( 2 a 2 - 1 )[ 2a 5 (5a 2 - 4a 2 +D “ “ 2 ! 

1 10( 2c2 2 - 1) ( 10o?2 ~ 1 2 c2 2 + 3) a a - ( 150q 2 - 260 a\ + 141q 2 - 24) 

60 ( 1 - a 2 ) (5a 2 - 2 ) ( 2 a 2 - 1 ) ( 1 - < 25 ) 


A 1_ 

° 5 ~ 60 


(5 a? - 2) (lOo?! ~ 12 q ? 2 + 3) 

a 5 (l - ce 5) (ce 2 - a 5 ) [2(5q ! 2 - 4g ? 2 + i)a 5 “ <^2] 


/ 
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12 . We still have to determine the coefficients /3 51 , (3 52 , ( 3 52 , Pu of our fifth- 
order formula. From equation (B) we find /3 51 . The fourth equations (Ai) 
and ( A2) together with equation (V,t) , determine the coefficients /3 52 , jS 53 , 
@M' 11 can b e shown that equation (V, t) is then also satisfied for any 
value of a 2 and . 

This concludes the computation of our RUNGE-KUTTA coefficients, since 
c q and c q can be determined from (I, 1) or (1, 1) , respectively, and the 

coefficients (3 (k = f , 2, 3, 4, 5) can be obtained from the standard equations 

K O 

K-l 

£ 0 . = a (K = 1,2, 3, 4, 5) (23) 

K A K 


13. Our RUNGE-KUTTA coefficients contain two arbitrary parameters, <a 2 and 
a 5 . We shall show that we can reduce the truncation error of our fourth- 
order formula by a proper choice of the parameter a 2 . From Table I, it 
follows that the leading term (the fifth -order term) of the truncation error 
consists of nine sub-terms. These sub-terms are certain expressions 
built up by the partial derivatives of the right-hand sides of the differential 
equation (1) . These sub-terms are multiplied by certain numerical fac- 
tors Tj through T 9 . For these factors, we find from equations V, f) 
through (V, 9) of Table I, replacing c by c and the upper limit 5 of the 
K-sums by 4 , 

T i = c 4 / 5 43^32 P 21 ~ 7377 \ 


T 2 - 


T 9 = 
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Naturally, it is desirable to find RUNGE-KUTTA formulas with small 
numerical factors Tj through T 9 to make the leading term of the local trun- 
cation error small. 

14. Since we can express all RUNGE-KUTTA coefficients that enter the right 

hand sides of (24) by a 2 alone, the factors T u . . . , T 9 can finally be written 
as functions of oi 2 - The computation results in the following values for these 

factors: 



2 

15. We see from (25) that all factors T 1? . . . , Tg would become zero for a 2 - g • 

Because of (20) , this value o; 2 leads to a 3 = 1. It can, however, be shown 
easily that a 3 = a 4 = 1 leads to contradictions in the equations (II, 1) , (III, 2) , 
(IV, 4) and (II, f) , (III, S) , (IV, %, (V,9). Therefore, we have to exclude 
2 

the value a 2 = ~ • 


Another interesting choice of ce 2 would result from 
10 q! 2 — 12 a? + 3 = 0 

If (26) would hold, all error factors in the second group of (25) would 
become zero. 


Because of (26) , it would follow from (21) , 


It can be shown easily that for the values ct 2 resulting from (26) and for 
c 4 = o, equation (B) cannot be satisfied. 
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Therefore, we also have to discard the values a 2 = (6 ±\I~6) 

resulting from ( 26) . 

However, by choosing for a 2 a value close to one of the above values, say 
close to 0. 355, we can expect that at least the error factors in the second 
group of ( 25) will become small. 


( 0,845 
^ 0. 355 


We shall consider two choices for a 2 that are reasonably close to 0.355 and 

1 

lead to relatively simple RUNGE-KUTTA coefficients, namely a 2 - 77 and 

3 ^ 

a 2 = ~ • The choice of a 5 in our formulas remains arbitrary. 

8 

16. Choosing o> 2 = ^ leads to the RUNGE-KUTTA coefficients of Table II. 


TABLE II. COEFFICIENTS FOR RK4( 5) , FORMULA 1 



Subtracting the last two columns of Table II from one another, one finds 
as approximation for the leading term of the local truncation error of our 
fourth -order formula 


12 





































TE 


f 

150 


o 


3 _ 16 1 . 

100 fz + 75 f 3 + 20 * 4 




(28) 


We also list the error factors through T 9 of our formula 1: 


1 


480 


4320 ’ 


T R = - 


480 ’ 


Te 4320 ’ 


T 9 " 17280 


T, =- 


T 7 = 


1440 


1 


5760 


T a - - 


4 4320 * 


2880 


>(29) 


17. For our second choice, — , we find the RUNGE-KUTTA coefficients 
of Table III. 


TABLE III. COEFFICIENTS FOR RK4( 5) , FORMULA 2 
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For the leading term of the local truncation error we obtain from Table III 
the approximation 


TE 


i 

360 o + 


128 f 2187 

4275 t2 + 75240 3 


S5* “ i f3 ) h (30) 


This formula 2 has the following error factors: 


rp _ ± rp __ x — hp _ 

1 _ 780’ 2 12480 ’ , 3 


1 


4160 ’ 


Ta = 


12480 


T R = - 


780 ’ 


12480 ’ 


T 7 - - 


1 


16640 


(31) 


T8 “ 8320 


Ta = - 


1 


49920 


We notice that the error factors (31) of our second formula are somewhat 
smaller than the corresponding terms (29) of our first formula. 

18. We should like to mention that another RK4(5) -formula was derived by 
D. SARAFYAN ( [3] , p. 4) . His fourth-order formula is based upon only 
four (instead of five) evaluations of the differential equations. Therefore, 
his fourth-order formula has larger error terms than our formulas. 


Since SARAFYAN's formula is published in an internal Technical Report 
and therefore is not easily accessible, we present SARAFYAN f s formula 
as Table IV. 

From Table IV it follows for the leading term of the local truncation error 


TE = (k + t f « + k - fk - S f 0 h 


(32) 


The error factors T 1 through T 9 for SARAFYAN’s formula read as follows 


i 


1 


Tr = 


Ta - 


■ 120 ’ 

1 

120 ’ 

1 

2880 


480 

1 


Te "480’ 


240 

1 


T, = - 


720 


T? 960 


Ts 480 ’ 


1 


(33) 
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TABLE IV. COEFFICIENTS FOR SARAFYAN’S RK4( 5) -FORMULA 



If we compare (33) with (31) , we notice that in (31) the factors T^ and T 5 
2 

are only — of the corresponding factors of (33) ; the other factors of (31) 

1 O 

are even smaller compared with the corresponding factors of (33) . 

For our RK4( 5) -formula number 1, the error factors are ~ (or better) of 
the corresponding error factors of SARAFYAN. 

Because of the smaller error factors we may expect that our RK4( 5) 
formulas 1 and 2 operate somewhat more economically than SARAFYAN’s 
formula. Numerical experiments on an electronic computer have confirmed 
these expectations. 


19. It is interesting to compare the error factors of SARAFYAN’s formula 
with those of KUTTA T s ( [4] , p. 443) standard fourth-order formula of 
Table V. 
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TABLE V. COEFFICIENTS FOR KUTTA’S RK 4- FOR MU LA 



The computation shows that KUTTA r s error factors are identical with the 
error factors (33) of SARAFYAN, except for T 7 , which in the case of 

KUTTA’s formulas has to be replaced by ■ . 

Therefore, SARAFYAN' s formula, in general, will not reduce substantially 
the number of integration steps required by KUTTA’s formula. However, 
it will reduce the computer time, since it requires only six evaluations per 
step compared to seven in KUTTA's formula, if the latter one is applied 
with the standard stepsize control procedure ( re computation of two steps 
as one step with double stepsize) . 

20. We might mention that we have presented general RUNGE-KUTTA trans- 
formation formulas with stepsize control in two earlier papers [5] , [6] . 

In the special case of a fourth-order formula these general formulas do not 
require any differentiation and can also be written in the form of classical 
RUNGE-KUTTA formulas. 

SECTION II. THIRD-ORDER FORMULAS 

21. Since a fourth-order RUNGE-KUTTA formula requires four evaluations 
(per step) of the differential equations, one might expect a pair of RUNGE- 
KUTTA formulas RK3(4) to require four evaluations also. 
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However, it can be shown easily that four evaluations are not sufficient to 
define the pair RK3(4) . The equations of condition lead to contradictions 
in the case of four evaluations. We therefore allow for five evaluations 
per step. 

However, it is possible to choose the fifth evaluation in such a way that this 
evaluation can be taken over as the first evaluation for the next step. There- 
by the number of evaluations per step again will be reduced to four, except 
for the very first step, when the integration is started. 

Since the derivation of the RK3(4) -formulas is very similar to the deriva- 
tion of the RK4( 5) -formulas of Section I, we may omit some details and 
present the main results only. 

22. The equations of condition, as they hold for a fourth-order formula, are 
listed in Table VI. 

TABLE VI. EQUATIONS OF CONDITION FOR FOURTH-ORDER FORMULA 
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TABLE VI. (Concluded) 


<IV, §) E 4 c k « k P k1 - | = 

K= 2 


A 1 \n’ A o 

(IV, 4) - • c a 3 - 


6 


K= i 


4a l 

K~K 24 


A similar table for a third -order formula is obtained from Table VI by 
omitting the fourth -order equations of condition and replacing, in the 
remaining equations, c by c and the upper limit 4 of the k - sums by 3. 

K K 

We denote the remaining four equations for the third -order RUNGE-KUTTA 
formula by (1,1), (11,1), (III, 1), and (III, 2). 

23. Again assuming that the conditions (5) hold and that d 2 , 03 are different 
from one another and from 0 and 1, we find from (II, 1) and (III, 2) : 


c 2 


3cr 3 - 2 
a 2 (o' 3 - d 2 ) 


c 3 - 


3d 2 - 2 
a 3 (a 2 - a 3 ) 


and from (II, 1) , (III, 2), (IV, 4): 


A 

C 2 - 


1 

12 


2a. 3 - 1 

G' 2 (0'3 " «2) ( 1 “ « 2 ) 


A 

C 3 = 


12 


2a 2 - 1 

a 3 (a 2 - d 3 ) ( 1 - a 3 ) 


A 

c 4 


1 6 q ! 2 d 3 - 4 (« 2 + d 3 ) + 3 

12 * (1 - a 2 )(l - d 3 ) 


(34) 




> (35) 


/ 


24. To bring the remaining equations of condition into a form that can be 

handled more easily, we make, similar to Section I, the further assump- 
tions 
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(B) C 2 /3 21 + C 3 /3 3! + 04/^4} = 0 

(D) /3 40 = c 0 , /3 41 = C} - 0 , /3 42 = c 2 , ^43 = c 3 

The assumption (6) is required if the fifth evaluation is to be taken over 
as the first evaluation for the next step. 


A A 

From the assumptions (A) , (B) , (D) , it follows immediately that the 
remaining equations of condition reduce to 


A A 


A 


(IV, 2) 

C3&32 a 2 

+ c 4 • 

3 


( £21 «i = 

1 2 

2“ 2 


(A) ' 

1 

(3zi a l + 

32^2 = 

f 2 

2« 3 

A 

A 

A 

0 

(B) 

c 2^2i + 

c 3 31 = 


12 


( 36 ) 


A 

25 . The first equation (A) expresses f$ 2 i a i by a 2- From equation (IV, 2 ) we 
obtain /3 32 as a function of a 2 and <23, since £ 3 and (^4 are given as functions 
of ce 2 and a 3 by ( 35 ) . 


The second equation (A) can then serve to find as function of a 2 and 

Qf 3. 

A 

The equation (B) , finally, represents a restrictive condition that must hold 
between 01 2 and a 3 to make the equations of condition compatible. The 
computation results in the following restrictive condition: 
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1 

2 


3q? 2 + i 


(37) 


«3 


3 « 2 


Eliminating a 3 from the coefficients of our third-order formula, these 
coefficients become functions of a 2 only. We find the following expressions 
for them: 


1 

1 

12 a? - 15a 2 + 

4 

° 2 " 6 

~af 

6 a 2 — 6 a 2 ^ 


2 

1 

(3a? - 2 ) (3«| - 3a 2 

+ l ) 2 

C 3 ~ 3 " 

JT 

Q ! 2 

0 

60^2 — 6 a? + 

1 

021«1 = 

1 2 

2“ 2 



031 a l = 

1 2 

( 3«2 “ 1 ) ( 3«2 “ 2 ) 


8«2 

(3aj> — 3 q ?2 + 1) ^ 



032 


1 6q?2 — 3q?2 + 1 


\ 


(38) 


/ 


26. We now consider the error factors for our third-order formula and try to 
make these error factors small by a proper choice of a 2 . 

From Table V we find the following four error factors: 

T l = c 3 /3 32 P2i “ 7 ^ 

T 2 = 7y( c 2 P 22 + C 3 P 32) 

T 3 = C 2 Ql 2^21 + c 3 a 3P31 

T 4 = ( c 2 a 2 + c&b “ 


24 

i 

8 


(39) 
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or, if we insert (38) into (39) : 


J_ (2 a 9 - 1) (3a? - i) 
24 3 g?2 *$<y,2 + i 


T _ * . ^t 2 ^ ~ 11 

Ts 8 1 3a2 - 3a 2 + i 

T _ i '(«; - 1) < 2a 2 - l) 2 
^ 8 3a 2 “ 2^2 + 1 


J_ (3g 9 - 1) (2g 9 - 1) 
24 3a 2 — 3a2 + i 


l (40) 


1 (a? ~ 1) ( 2 <^:> ~ 1) 2 

24 3a 2 — 3 a 2 + i 


From the second equation (40) , it follows that we can make T 2 = 0 by 
choosing 

= -32^1 (41) 

1 3a 2 

27. All four error factors would become zero for a 2 = j • This means for this 

choice of a 2 , our third-order formula would actually become a fourth-order 
formula. Our pair RK3(4) of RUNGE-KUTTA formulas of the third- and of 
the fourth -order would degenerate into one fourth -order formula. 

However, by choosing for a 2 a value close to — , we obtain RK3(4) -formu- 
las with small error factors (40) . We give two examples for such formulas 
with small error factors. 


28. Choosing a 2 = “ , we obtain the RUNGE-KUTTA coefficients of Table VII. 

9 

For the leading term of the local truncation error we find from Table VII 


TE - 


27 


. _245_ f \ . 

U + . f 3 12 f4 ) h 


288 "o 416 ' 1872 ^ 

For the four error factors of the formula of Table VII we obtain 

5 5 


(42) 


Ti - 


1 

168 ’ 


t 2 = 0 , t 3 


1512 


4536 


(43) 
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For the formula of Table VIII we obtain 


/ 4 75 P 5776 , 1 , \ , 

TE = ^735 f o “ 3626 ^ + 81585^ " 18 4 j 

and the four error factors 

Tl = 228 s T 2 - 0 , T 3 - — , t 4 2565 


( 44 ) 


(45) 


30. For comparison, we list KUTTA’s third-order formula ([4], p. 440) in 
Table IX. His formula reads: 


TABLE IX. COEFFICIENTS FOR KUTTA’S RK3-FORMULA 



It has the following error factors: 

T i = - L, t 2 = 0, T 3 = 

Comparing (43) and (45) with (46), we 

is only — and the largest value of (45) 

(46). 7 

Therefore, our formulas of Table VII and Table VIII can be expected to 
be more economical than KUTTA’s RK3 -formula. 

Since KUTTA’s formula does not provide for a stepsize control, it requires 
five evaluations per step if operated with the standard stepsize procedure 
(recomputation of two steps as one step with double stepsize) . 


T 4 = 0 (46) 

see that the largest value of (43) 
o 

only — of the largest value of 
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SECTION III. SECOND-ORDER FORMULAS 


31. Allowing for four evaluations per step for an RK2(3) -formula, we list 
in Table X the equations of condition for a third-order formula. 

TABLE X. EQUATIONS OF CONDITION FOR THIRD-ORDER FORMULA 


A 

(I,D 


E 

K= 0 


c -1=0 

K 


A 

(n, i) 


E 

K=1 


A 

c a 
K K 


2 = 0 


( 111 , 1 ) 


z 

K= 2 


c P , 

K K 1 6 


- - = 0 


(III, 2) 


1 A 2 

2 ^ C K a K 6 

K - 1 


- - = 0 


32. We want to use the fourth evaluation of the differential equations as 
the first evaluation for the next step. This requires the conditions 

a 3 - 1 , £ 30 " c 0 > 031 = C 1 » 032 ~ c 2 ( 47 ) 

where the c r s are the weight factors of the second-order formula, which is 

v 

obtained from the first two equations of Table IX, replacing c by c^ and 

the upper index 3 of k- sums by 2. We denote those equations for the 
second-order formula by (I, 1) and (II, 1) . 

Furthermore, we assume 

c t = 0 (48) 
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and 


(A) 


21 


1 2 

T«2 


P 31 “ o 


33. Because of (A), equations (III,i) and (III, 2) become identical. Therefore, 
we can omit equation (III, 1) from Table X. 

A A A , 

From equations (II, 1) and (III, 2) we obtain, because of c t = 0 and a 3 = 1, 

1 


A 1 

C 2 - “• 


A 

c 3 


6 a 2 ( I " ^ 2 ) 

1 . 2 - 3q> 2 
6 1 - a? 


(49) 


34. Let us now consider the error factors of our second -order formula. From 
the last two equations of Table X, we obtain for the two error factors the 
expressions 

1 

T 1 = c 2 P2i a l “ 6 

T 2 = ^ ( C l a l + C 2 Q! 2 ) - “ 

Because of the first equation (A) , we 

rp f 2 ^ rp 

Ti - i C 2“2 - q = T 2 

or , 

T 2 - Tj = < 52 ) 

From (52), it follows that T 2 * T t , assuming c t * 0 . 


(50) 


can write for the first equation ( 50) 


1 2 

7T c l«l 


(51) 
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For the following, we might assume 


t 2 =- t i (53) 

By a proper choice of and a 2 we try to make T* and T 2 = - T 1 sufficiently 
small. 


From (51), (52), and (53) we find 

rp 1 2 f 1 9 

Ti = g c 2 a 2 - - = - -c ia t 
or 

2 

c i<*i + 2 c 2 q;2 = 


(54) 


Equation ( 54) and equation (II, 2) represent a system of two linear equa- 
tions for Cj and c 2 . 


The system has the solution 

i 2 “ 3 q! 2 
c 1 — • * 

3 Q! 1 (q' 1 - 2of 2 ) 

i . 3ogi - 4 

° 2 6 oizioii - 2a 2 ) 


(55) 


Introducing the expression (55) for c 2 into (51) yields 


t - 3q 2 - 2 

1 12 1 Q!j - 2d!2 


(56) 


From equation ( 56) , it follows that we could make Tj — T 2 — 0 by choosing 

2 

ol 2 = g- . However, as in the case of our third-order formula RK3(4) , our 

pair of RUNGE-KUTTA formulas RK2(3) would then degenerate into a 
single third -order formula. 

2 

However, by choosing a 2 close to ~, we might obtain a suitable pair 
of RK2(3) -formulas with small error factors T t and T 2 . 



Choosing a 2 = and a t = j , we find the RUNGE-KUTTA coefficients 
of Table XI as follows . 



From Table XI we find for the leading term of the local truncation error 
the approximation 


TE = 


23 J_ 

1782*° + 33 1 


350 , 

fo + 

11583 1 



h 


The error factors for our RK2{3) -formula would read 


Ti 


1 

2112 ’ 


t 2 


1 

2112 


(57) 


(58) 


36. It is well-known that second -order RUNGE-KUTTA formulas RK2 can be 
obtained by two evaluations only. It is also possible to provide a steps ize 
control procedure by a third evaluation. We list in Table XII an example 
for such a RK2(3) -formula. 

For the truncation error we find from Table XII the approximation 


TE = 




(59) 
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TABLE XII. COEFFICIENTS FOR RK2(3) , BASED ON THREE EVALUATIONS 



The error factors for the formula of Table XII would read 



It is interesting to note that KUTTA’s third-order formula (Table IX) can 

also be operated as a second-order formula RK2(3) with steps ize control. 

We have to change only the weight factors of the formula, considering the 

c -values of Table IX as c -values of the RK2(3) -formula and using c 0 = D 
K K 

Ci = 1 as its c -values. 

1 K 

Comparing (58) with (60) , we notice that the error factors of our formula 
of Table XI are only of the larger error factor (60) . Therefore, we 

O Ow 

may again expect our formula of Table XI to be more economical than the 
formula of Table XII, since the formula of Table XI also requires only three 
evaluations per step (except for the very first integration step) . 


SECTION IV. FIRST-ORDER FORMULA 


37. If we base our RUNGE-KUTTA formula RKi ( 2) upon three evaluations per 
step, we obtain for the second -order formula the following equations of 
condition. 
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TABLE XIII. EQUATIONS OF CONDITION FOR SECOND-ORDER FORMULA 


A 

(i,D 


z 2 

K=0 


A 

C 

K 


1 = 0 


A 

(ii, i) 


z 

K — 1 


A 

c a 

K K 


t 

2 


= 0 


Since we intend to use the third evaluation again as the first evaluation for 
the next step, we require 


OL 2 — 1 , fi 20 - c 0 > @21 — C 1 


(61) 


where c 0 and c 1 are the weight factors of the first-order formula, which is 
obtained from the first equation of Table XIII, replacing c^ by c^ and the 

upper index 2 of the k- sum by 1. Let us denote this equation for the first- 
order formula by (1, 1) . 

A 

Obviously, there is only one error factor T 1? obtained from (II, i) as: 


1 

Ti - Citti - - 


(62) 


We have to make this error factor small to obtain an efficient RK1(2) - 
formula. However, we should not make Tj zero, since our pair of formulas 
RK1 ( 2) would degenerate for T t = 0 into one second-order formula. We 
choose as coefficients the values of Table XIV: 



which clearly satisfy equations (I, 1) , (1, 1 ), an d (II, i ). 
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For the approximate truncation error we find 


TE = ^(fo-f 2 )h, (63) 

and the error factor Tj becomes 



38. EULEK-CAUCHY T s method can be considered as a first-order RUNGE- 
KUTTA formula, requiring only one evaluation per step. The method can 
even be operated as a RUNGE-KUTTA formula RK1(2) without additional 
evaluations by adding a second evaluation that can again be taken over as the 
first evaluation for the next step. 

The resulting pair of formulas is shown in Table XV. 

TABLE XV. EULER-CAUCHY’S METHOD AS RK1(2) 



The second-order formula RK2 of Table XV is obviously the so-called 
modified EULER-CAUCHY method. 

Table XV yields as approximate truncation error 

TE = |(f 0 -fi)h (65) 

The error factor for the RKi(2) -formula of Table XV becomes 
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Comparison of (64) and (66) clearly shows that our formula of Table XIV 
will be more efficient than the formula of Table XV. Even the one additional 
evaluation per step required by our formula of Table XIV, in general, will 
not outweigh the advantage of the much smaller error factor. 

We shall present a numerical example in Section V. 


SECTION V. A NUMERICAL EXAMPLE 
(ORDINARY DIFFERENTIAL EQUATIONS) 


39. In this section we present the numerical results of our new formulas for 

the same example that we considered in our earlier NASA report ([!],. p. 30) 

y' - - 2xy • logz , z ? = 2xz 9 logy \ 

Initial values: x 0 = 0 , y 0 = e , z 0 = 1 1 (67) 

, cos(x 2 ) sin(x 2 ) I 

Exact solution: y = e , z = e / 

For reasons of comparison, we also include the results for some formulas 
of other authors. Our results are presented in Table XVI, the results of 
our new formulas being marked by an asterisk (*) . 

The numerical integration of problem (67) was executed on an IBM-7094 
computer. All methods listed in Table XVI were run in single precision 
(eight decimal digits) and with the same tolerance (10 8 ) for the local 
truncation error. Since the first-order formulas require a very small 
stepsize, we ran these formulas only up to x = 5; for the higher -order 
formulas the integration was performed up to x - 25. 

40. Table XVI shows that the small error factors of our new formulas (*) 
make themselves felt in a considerable reduction of the number of integra- 
tion steps and, therefore, of the computer time required for our problem. 

The greatest reduction gained was for our first- and second-order formulas. 

1 1 

For these formulas we required only about — and about — of the computer 

time of the corresponding conventional formulas. In the case of the third- 

1 

order formula, we could cut the computer time to less than — , compared 
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TABLE XVI. COMPARISON OF THE VARIOUS METHODS FOR EXAMPLE (67) 

co 

to 


Results for x = x (Tolerance: 10" 8 ) 
max 


Method 

Order 

of 

Method 

Number 

of 

Evaluations 
Per Step 

X 

max 

Number 

of 

Steps 

Total 

Number 

of 

Evaluations 

Running- 

Time 

(min) 

on 

IBM-7094 

Accumulated Er 
Ay 

rors in y and z 
A z 

EULER-CAUCHY (Table XV) 
RK1(2) (Table XIV) 

1 

1 

1 

2 

5 

5 

269 956 
16 871 

269 956 
33 742 

2. 90 
0.32 

+0. 3018 • 10" 2 
+0. 1926 • 10“ 3 

-0. 2945 • 10‘ 3 
-0. 1543 • 10 -4 

RK2(3) (Table XII) 
RK2( 3) (Table XI) 

2 

2 

3 

3 

25 

25 

243 510 
37 493 

730 530 
112 479 

6.43 

1.13 

+0. 1458 • 10' 4 
-0. 1874 • 10” 4 

-0. 1781 • 10“ 3 
-0. 8330 * 10“ 5 

KUTTA ( Table IX) 
RK3( 4) (Table VII) 
RK3( 4) (Table VIII) 

3 

3 

3 

5 

4 

4 

25 

25 

25 

41 862 
23 225 
22 054 

209 310 
92 900 
88 216 

1.91 

0. 90 
0. 85 

+0.2190 • 10" 5 
-0.2611 • 10" 5 
-0.2578 • 10" 5 

-0. 4664 • 10~ 5 
+0. 1639 • 10" 4 
+0. 1474 • 10 -4 

KUTTA ( Table V) 
SARAFYAN (Table IV) 
RK4( 5) ( Table II) 

RK4( 5) ( Table III) 

4 

4 

4 

4 

7 

6 

6 

6 

25 

25 

25 

25 

16 010 
14 746 
11 059 
9 947 

112 070 
88 476 
66 354 
59 682 

0.95 
0. 81 
0.68 
0.58 

+ 0. 1881 • 10" 5 
-0. 1546 • 10 -5 
+0. 1222 • 10 -5 
+0.2041 • 10“ 5 

+0.2207 • 10" 4 
-0.2086 • 10"° 
+0.2015 • 10“ 4 
+0.2512 • 10" 4 


(.*) 






with KUTTA’s formula. For fourth-order formulas our gains are more 
modest. Out second fourth-order formula (Table III), however, still 

3 

requires only about 7 - of the computer time of KUTTA’s formula and about 

Q 

— of the time of SARAFYAN’s formula. 

4 

The accuracy of our new formulas is about the same as that of the conven- 
tional formulas, except for the case of the first-order formulas. In this 
case we gain one more decimal digit with our new formula. Because the 

number of steps required by our new formula is only about — of the 

number of steps of EULER-CAUCHY’s formula, we do not accumulate 
round-off errors as heavily as the latter formula. 

Compared with the conventional RUNGE-KUTTA formulas our new formulas 
offer results of the same accuracy in a fraction of the computer time. 
Therefore, our new formulas might be of interest for the numerical inte- 
gration of ordinary differential equations. 

Moreover, we shall show in the following sections of this paper that they 
can also be applied successfully to certain partial differential equations of 
the parabolic type. 


SECTION VI. APPLICATION TO HEAT TRANSFER PROBLEMS 


41. Let us consider the one -dimensional heat transfer problem 


9 u _ 9 2 u 

9 t ” 9? 

Initial Condition: t 


Boundary Conditions: 

- 


x 

x 


= 0 : u(x, 0 ) 

= 0 : u( 0 , t) 

= 1 : u ( 1 , t) 


u 0 (x) 

b 0 ( t) 
b t (t) 


( 68 ) 


The first equation of ( 68 ) represents the simplest partial differential equa- 
tion of the parabolic type. 
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Problem (68) requires finding a quantity u (usually the temperature) as 
function of the space variable x and the time variable to 

The existence of a solution of the problem (68) can be shown, under proper 
assumptions, by Fourier series methods. See, for example, the textbook 
by CARSLAW and JAEGER ( [7] , pp. 76-88) „ A quite elementary proof of 
the uniqueness of the solution of (68) is given in a textbook by BIEBERBACH 
( [8] , pp. 352-353) . 

42. Replacing both derivatives of the first equation (68) by finite differences, 
one obtains the well-known difference equation 


u. . . - u. . 

Ej + 1 j 


k 


u 


i± l > J 


2 u. + 

"hT 


u 


Rid 


(69) 


with k = At and h = Ax . 

Equation (69) is widely used for obtaining numerical solutions of problem 
(68) , since this is the simplest explicit approach to the problem. 

In applying (69) one has to pay attention to the fact that the mesh sizes h 
and k have to satisfy the well-known stability condition 



The condition ( 70) represents a certain restriction for the time mesh k 
if the space mesh h is given. 

To preserve a reasonable accuracy of the results one has to apply a small 
space mesh h. In some problems the time mesh k, resulting from (70) , 
may then become prohibitively small. 

Furthermore, one has to consider that the mesh size k, resulting from 
( 70) , will guarantee the stability of the results, but no estimate of the 
local truncation error can be obtained from ( 70) . 

43. There exist other difference methods without stability restrictions, for 
example, the explicit method of DU FORT-FRANKEL [9] : 
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, (71) 


u. 


L 3+1 


- u 


kiri 


i+k j 


- (u. 


kill 


u 


k j" 1 


) + u 


i-k j 


2k 


IT 


or the implicit method of CRANK-NICOLSON [10]: 


u. . - u. . 

L.1+1 LJ 

k 


( u , " 2u. + u 

1 J i+1 , j+1 L.1 + 1 *~k3 + 1 

2 I h 7 ” 


u , . - 2u. . + u . 
i+l,3 hi iliii 
+ Y? 


72) 


However, other difficulties occur in the application of these two methods. 
The method of DU FORT-FRANKEL is a three-level method (time levels 
j ^ j+i) requiring a special starting procedure and additional consid- 
erations for a stepsize change of the time step. In the case of the implicit 
method of CRANK-NICOLSON, one has to solve a system of linear equa- 
tions of the triple -diagonal form for each time step. 

Although these difficulties might not be considered too serious, these 
methods again do not furnish an estimate of the local truncation error. 

Since they have no stability restrictions, an increased danger exists in 
that the integration could be performed with too large a time step. If this 
is the case, one would lose any accuracy, and, after a certain number of 
time steps, the computed u— values would bear little resemblence to the 
solution of the problem. 

Equations (72) and (69) have identical left-hand sides; they approximate the 

time derivative by the same first-order accurate, finite difference. 

8 t 

Therefore, the time integration of the CRANK-NICOLSON method cannot 
be more accurate than in the case of the explicit method (69) . 

44. If one wants to take full advantage of the accuracy that modern electronic 
computers are capable of, a more accurate approximation of the partial 
derivatives in the heat transfer equation first equation (68) and a 
current stepsize control for the integration become absolute necessities. 

A convenient way to achieve these objectives consists of converting the 
heat transfer equation into a system of ordinary differential equations 
by replacing the space derivative only by a finite difference . The resulting 
system of ordinary differential equations can then be integrated by RUNGE 
KUTTA methods. 
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This is not a new procedure; it was suggested long ago by D. R. HARTREE 
and collaborators ( [11] and [12], p. 254) . They suggest replacing the first 
equation (68) by the system 


du.(r) 
d t 


u i+ i (T) 


2u(t)+ u (t) 


6 2 u .. ( t ) 


(73) 


with r = ty . 

hr 

0 U 

Naturally, one can apply a better approximation to by using higher- 

order central differences, resulting in 

d U (t) 

— t = 5 2 u. (r) - — 6 4 u.(t) + — 5 6 u. ( r) (74) 

dr i 12 1 90 1 

Formula (74) is suggested in GOODWIN*s book ( [ 13] , p. 113) . 

In equations (73) and (74) the authors used the standard fourth-order 
RUNGE-KUTTA formula (Table V) for the numerical integration. 

GOODWIN ( [13] , pp. 114-115) and FOX ( [14] , pp. 240-241) point out that 
the fourth -order RUNGE-KUTTA method offers few advantages compared 
with the explicit difference formula (69) , since the stability restriction 
of the RUNGE-KUTTA method is only slightly better (r < 0. 7 instead of 
r < 0. 5) . On the other hand, the RUNGE-KUTTA method requires con- 
siderably more computer time per time step than formula (69) . There- 
fore, it might seem questionable whether RUNGE-KUTTA methods are 
really worthwhile for application to problems such as (68) . 

45. We feel that the explicit difference method (69) certainly has its merits 
when applied to the computationally simple problem (68) and when requir- 
ing a moderate accuracy only. 

However, the new low-order RUNGE-KUTTA formulas withstepsize control 
that are presented in this paper show the application of RUNGE-KUTTA 
formulas to parabolic partial differential equations in a different light, 
especially in connection with the use of highly accurate electronic computers. 

For our numerical computations we used an eight -decimal digit computer 
(IBM-7094) , and, when applying our new RUNGE-KUTTA formulas to heat 
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transfer problems, we required that the local truncation error became 
negligible. This means we required that the local truncation error 
not contribute to the eight leading digits that the computer carries. 

By doing so, we were able, when using higher-order central differences 
for the space derivatives, to obtain by means of our new RUNGE-KUTTA 
formulas, five to six good digits, even after a considerable number of 
integration steps. ( See the examples in Section VII. ) We found that these 
somewhat severe accuracy requirements could be well -satisfied by our 
low-order RUNGE-KUTTA methods and that we did not run into stability 
problems with these low-order methods since our accuracy requirements 
are much sharper than the stability requirements of our formulas. 

46. Comparing the explicit difference method (69) with EULER-CAUCHY’s 
method (Table XV) , we see that both methods are identical as far as the 
time integration is concerned. Since the method of Table XV yields a 
convenient stepsize control procedure we have substituted EULER-CAUCHY’s 
method of Table XV, for the explicit difference method (69) . 

47. Naturally, the explicit difference method and any RUNGE-KUTTA method 
can be run in a fraction of the computer time if we apply less severe 
accuracy restrictions, for example, if we apply the stability restriction 
only. However, such a relaxation of the accuracy restriction would result 
in a severe loss of accuracy. 

48. The RUNGE-KUTTA methods are not restricted to the simple problem (68) . 
They can be applied to more involved one -dimensional heat transfer prob- 
lems, as outlined in the next section. Multi -dimensional heat transfer 
problems can also be made amenable to RUNGE-KUTTA methods when 
replacing all space derivatives by proper finite differences. 


SECTION VII. TWO NUMERICAL EXAMPLES 
(HEAT TRANSFER PROBLEMS) 


49. In this section we present numerical results of our new RUNGE-KUTTA 
formulas for two one -dimensional heat transfer problems. In both prob- 
lems the exact solution is available so that the errors of our formulas can 
be stated. 
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First problem: 


9 u 


v 2 


-u 


9 t 4 2 + x 

Initial Condition: 

Boundary Conditions: 


9 2 u 
9 x' 2 


t - 0: u = 2^1- log( 2 - x 2 )] 


( x - 0 : 


9 u 
9 x 


x = 1 : u - 2 + log( 1 4- t) 

Exact Solution: u = 2 + log( 1 + t) - 2 log( 2 - x 2 ) 


( 75 ) 


9 ^ u 

Using higher-order central differences when replacing - — jr- by finite 

O X 

differences, the partial differential equation of our problem is converted 
into the following system of ordinary differential equations: 


du^T) 

d r 


i 

4 


2 + x 


l 


'V T) i . i 

e J 6 2 u.(t) - — <5 4 u.(t) 

+ h 6 \ (T) _+ - • ' ) 


(76) 


The introduction of the fourth-, sixth-, . . . order differences in ( 76) 
makes necessary some extrapolation in the vicinity of the boundary line 
x - 1 . When including the fourth-order differences only, we have assumed 
that the sixth-order differences are constant in tie vicinity of x - i. 
Correspondingly, the eighth-order differences are assumed to be constant 
if sixth-order differences are still carried in ( 76) . For x = 0 no extrap- 
olation is needed since the problem is symmetric with respect to x = 0. 


For our computation we divided the x-interval < 0, 1 > into 16 equal parts 

(h = ~ ) . The numerical integration of ( 76) with respect to r = -77 
16 h 

was performed with a tolerance of 10 -8 on an eight-decimal digit computer. 

The steps ize control procedure was applied for x = 0 only. 


Table XVII shows the results for t - 100, obtained with some of our new 
formulas. Under the heading "Method” of Table XVII we have also indicated 
the order of the highest central difference 6 2 , , or <5 8 carried in ( 76) . 
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TABLE XVII. COMPARISON OF VARIOUS RUNGE-KUTTA METHODS 

FOR PROBLEM { 75) 


Results for t = 100 (Tolerance: 10“ 8 ; h = — 


“O 



Order 



Number 

Running Time 

X. 

Accumulated Error 


of 

of 

(min) 

(maximum) 

Method 

Method 

Steps 

on IBM-7094 

in u 

EULER-CAUCHY (Table XV) - 6 2 

1 

30 721 

1. 75 

0. 1408 • 10 -2 

RK1 ( 2) (Table XIV) - <5 2 

1 

1 924 

0.22 

0. 1452 • 10- 2 

RK2( 3) (Table XI) - S 2 

2 

822 

0.23 

0. 1425 * 10~ 2 

RK3 ( 4) (Table VIII) — 6 2 

3 

1 036 

0.33 

0. 1424 . 10 -2 

EULER-CAUCHY (Table XV) - <5 4 

1 

30 715 

2. 07 

-0. 3767 • 10" 4 

RK1 ( 2) (Table XIV) — <5 4 

1 

1 998 

0.28 

-0. 1991 • 10 -4 

RK2( 3) (Table XI) - <5 4 

2 

1 070 

0.36 

-0.1961 • 10 -4 

RK3( 4) (Table VIII) — 6 4 
— 

3 

1 305 

0. 50 

-0. 2086 • 10~ 4 

EULER-CAUCHY (Table XV) -6 e 

1 

30 712 

2. 33 

-0.2068 • 10 -4 

RK1 (2) (Table XIV) - 6 e 

1 

2 072 

0. 39 

-0. 1848 • 10~ 5 

RK2( 3) (Table XI) - S 6 

2 

1 227 

0.45 

-0.3636 • 10~ 5 

RK3 ( 4) (Table VHI) — o 6 

3 

1 520 

0.61 

-0. 4947 • 10~ 5 


From Table XVII one recognizes immediately the gain of accuracy obtained 
by taking into account the fourth- and sixth-order central differences in 
( 76) . The results clearly suggest that at least the fourth-order central 
differences should be included in the computation. 

The table also shows that our new RUNGE-KUTTA methods give about the 
same or more accurate results in a fraction of the time required by EULER - 
CAUCHY f s method, which is identical with the explicit difference formula 
(69) . 

However, our third-order formula is already slower on the computer than 
our first- and second-order formulas. The trend continues: our fourth- 
order formula is again slower than our third-order formula, etc. 

This confirms the experience that other authors have had with fourth - 
order RUNGE-KUTTA formulas, which are too slow for heat transfer 
problems and cannot compete with our new lower-order RUNGE-KUTTA 
formulas . 
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50. We present another example: 


Second problem: 

9 u _ 9 2 u 

9 t ~ 9 x 2 


+ 2t 


+ u [ ( logu) 2 + log u 


1 ] 


Initial Condition: 


Boundary Conditions: 


Exact Solution: 


t = 

x 

< 

X 

u = 


cosx 


0 : u = e 

= 0 : u = e 

- 1 : u = e 

cos(x+t 2 ) 

e 


cos(t 2 ) 
cos( 1+t 2 ) 




> ( 77 ) 


0 

In this example we have to replace — by finite differences also. Taking 

into account higher-order central differences, we replace the partial dif- 
ferential equation of our problem with the following system of ordinary 
differential equations: 


du (t) 
d r 



The numerical integration of ( 78) is performed in the same way as in our 
first problem. However, our new problem requires some extrapolating 
for x = 0 also; and the stepsize control procedure is now performed for 

i 
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Table XVIII shows our results for t = 5. These results are quite similar 
to those for our first problem. The computer time for our new RUNGE- 
KUTTA formulas is, however, here a considerably smaller fraction of the 
computer time of the explicit difference formula ( EULER-CAUCHY's 
formula) . 

TABLE XVIII. COMPARISON OF VARIOUS RUNGE-KUTTA METHODS 

FOR PROBLEM (77) 

Results for t = 5 (Tolerance: 10~ 8 ; h = — ) 

' 16 ' 




Order 


Running Time 

Accumulated Error 



of 


(min) 

(maximum) 


Method 

Method 


on IBM-7094 

in u 

EULER- 

-CAUCHY (Table XV) - o 2 

1 

235 354 

19.46 

0.6313 • 10“ 3 

RK1 ( 2) 

(Table XIV) — 6 2 

1 

14 737 

2. 47 

0. 6711 • 10 3 

RK2( 3) 

(Table XI) - 6 2 

2 

2 142 

0. 87 

0. 7065 • 10~ 3 

RK3( 4) 

(Table VIII) - 6 2 

3 

2 519 

1. 03 

0. 6745 • 10" 3 

EULER 

-CAUCHY (Table XV - 6 4 

1 

235 388 

22. 87 

-0.4032 • 10' 4 

RK1 ( 2) 

(Table XIV) - 6 4 

1 

14 906 

2. 93 

-0.3516 ■ 10~ 5 

RK2( 3) 

(Table XI) - 6 4 

2 

2 695 

1.08 

0.3994 • 10 -5 

RK3( 4) 

(Table VIII) - 6 4 

3 

3 334 

1. 91 

0. 3129 • 10 -5 

EULER 

-CAUCHY (Table XV) - 6 6 

1 

235 372 

25. 32 

-0.5391 • 10" 5 

RK1 ( 2) 

(Table XIV) - b 6 

1 

14 971 

3.25 

-0.4351 • 10~ J 

RK2( 3) 

(Table XI) - 6 8 

2 

2 812 

1. 50 

-0. 3427 • lO" 5 

RK3 ( 4) 

(Table VIII) - < 5 6 

3 

3 677 

2.38 

-0. 1132 • 10 5 


51. In both examples, our new low-order RUNGE-KUTTA formulas proved to 
be considerably faster and of equal or better accuracy than the conventional 
explicit difference formula (69) . This is not too suprising if one realizes 
that (69) is a formula of first-order accuracy only, and that its truncation 
error factor is 256 times as large as the error factor of our first-order 
formula RKi(2) of Table XIV. 


George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama 35812, April 15, 1969 
129-04-03-00-62 
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